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Quantum-critical strongly correlated electron systems are predicted to feature universal collision- 
dominated transport resembling that of viscous fluids. nmi However, investigation of these phenom¬ 
ena has been hampered by the lack of known macroscopic signatures of electron viscosity [&H9j. Here 
we identify vorticity as such a signature and link it with a readily verifiable striking macroscopic 
DC transport behavior. Produced by the viscous flow, vorticity can drive electric current against an 
applied field, resulting in a negative nonlocal voltage. We argue that the latter may play the same 
role for the viscous regime as zero electrical resistance does for superconductivity. Besides offering 
a diagnostic which distinguishes viscous transport from ohmic currents, the sign-changing electri¬ 
cal response affords a robust tool for directly measuring the viscosity-to-resistivity ratio. Strongly 
interacting electron-hole plasma in high-mobility graphene I10lll2| affords a unique link between 
quantum-critical electron transport and the wealth of fluid mechanics phenomena. 


Symmetries and respective conservation laws play a 
central role in developing our understanding of strongly 
interacting states of matter. This is the case, in particu¬ 
lar, for many systems of current interest, ranging from 
quantum-critical states in solids and ultracold atomic 
gases to quark-gluon plasmas PJ0, which share com¬ 
mon long-wavelength behavior originating from the fun¬ 
damental symmetries of space-time. The ensuing energy 
and momentum conservation laws take the central stage 
in these developments, defining hydrodynamics that re¬ 
veals the universal collective behavior. Powerful ap¬ 
proaches based on conformal field theory and AdS/CFT 
duality grant the well-established notions of fluid me¬ 
chanics, such as viscosity and vorticity, an entirely new 
dimension [13l [14] . 

Despite their prominence and new paradigmatic role, 
viscous flows in strongly correlated systems have so far 
lacked directly verifiable macroscopic transport signa¬ 
tures. Surprisingly, this has been the case even for con¬ 
densed matter systems where a wide variety of experi¬ 
mental techniques is available to probe collective behav¬ 
iors. Identifying a signature that would do to viscous 
flows what zero electrical resistance did to superconduc¬ 
tivity has remained an outstanding problem. The goal 
of this article is to point out that vorticity generated in 
viscous flows leads to a unique macroscopic transport be¬ 
havior that can serve as an unambiguous diagnostic of the 
viscous regime. Namely, we predict that vorticity of the 
shear flows generated by viscosity can result in a back- 
flow of electrical current that can run against the applied 
field , see Fig0 The resulting negative nonlocal voltage 
therefore provides a clear signature of the collective vis¬ 
cous behavior. Associated with it are characteristic sign¬ 
changing spatial patterns of electric potential (see Fig@ 
and Fig [2]) which can be used to directly image vortic¬ 
ity and shear flows with modern scanning capacitance 
microscopy techniques. [131 

The negative electrical response, which is illustrated 


in Fig jT] originates from basic properties of shear flows. 
We recall that the collective behavior of viscous systems 
results from momenta rapidly exchanged in carrier col¬ 
lisions while maintaining the net momentum conserved. 
Since momentum remains a conserved quantity collec¬ 
tively, it gives rise to a hydrodynamic momentum trans¬ 
port mode. Namely, momentum flows in space, diffusing 
transversely to the source-drain current flow and away 
from the nominal current path. A shear flow established 
as a result of this process generates vorticity and (for 
an incompressible fluid) a back flow in the direction re¬ 
verse to the applied field. Such a complex and manifestly 
non-potential flow pattern has a direct impact on the 
electrical response, producing a reverse electric field act¬ 
ing opposite to the field driving the source drain-current 
(see Figj2]). This results in a negative nonlocal resistance 
which persists even in the presence of fairly significant 
ohmic currents (see Fig{2]). 

Attempts to connect electron theory with fluid me¬ 
chanics have a long and interesting history, partially 
summarized in Refs. [8, 18, 19 . Early work on viscos¬ 
ity of Fermi liquids made connection with ultrasound 
damping. m Subsequently, Gurzhi introduced an elec¬ 
tronic analog of Poiseuille flow. m Related tempera¬ 
ture dependent phenomena in nonlinear transport were 
observed by deJong and Molenkamn.[l8] Recent devel¬ 
opments started with the theory of a hydrodynamic, 
collision-dominated quantum-critical regime advanced by 
Darnle and Sachdev.pQ Andreev, Kivelson, and Spivak 
argued that hydrodynamic contributions can dominate 
resistivity in systems with a large disorder correlation 
length. [7] Forcella et al. predicted that electron viscosity 
can impact electromagnetic field penetration in a dra¬ 
matic way. 0 Davison et al. linked electron viscosity 
to linear resistivity of the normal state of the copper 
oxides. [20] 

As a parallel development, recently there was a surge of 
interest in electron viscosity of graphene. 0 m n in I2z3 
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FIG. 1: Current streamlines and potential map for vis¬ 
cous and ohmic flows. White lines show current stream¬ 
lines, colors show electrical potential, arrows show the direc¬ 
tion of current. Panel a) presents the mechanism of a negative 
electrical response: Viscous shear flow generates vorticity and 
a back flow on the side of the main current path, which leads 
to charge buildup of the sign opposing the flow and results in 
a negative nonlocal voltage. Streamlines and electrical poten¬ 
tial are obtained from Eq. § and Eq.|6|. The resulting po¬ 
tential profile exhibits multiple sign changes and ±45° nodal 
lines, see Eq. 0- This provides directly measurable signatures 
of shear flows and vorticity. Panel b) shows that, in contrast, 
ohmic currents flow down the potential gradient, producing a 
nonlocal voltage in the flow direction. 


The quantum-critical behavior is predicted to be par¬ 
ticularly prominent in graphene. j 1 OUTS] Electron inter¬ 
actions in graphene are strengthened near charge neu¬ 
trality (CN) due to the lack of screening at low car¬ 
rier densities. ns ca As a result, carrier collisions are 
expected to dominate transport in pristine graphene in 
a wide range of temperatures and dopings. [251 Further¬ 
more, estimates of electronic viscosity near CN yield one 
of the lowest known viscosity-to-entropy ratios which ap¬ 
proaches the universal AdS/CFT bound. [5] 

Despite the general agreement that graphene holds the 
key to electron viscosity, experimental progress has been 
hampered by the lack of easily discernible signatures in 
macroscopic transport. Several striking effects have been 
predicted, such as vortex shedding in the preturbulent 
regime induced by strong current [B], as well as nonsta¬ 
tionary flow in a ‘viscometer’ comprised of an AC-driven 
Corbino disc. [9] These proposals, however, rely on fairly 
complex AC phenomena originating from high-frequency 
dynamics in the electron system. In each of these cases, 
as well as in those of Refs. ® uni, a model-dependent 
analysis was required to delineate the effects of viscosity 
from ‘extraneous’ contributions. In contrast, the nonlo¬ 
cal DC response considered here is a direct manifestation 
of the collective momentum transport mode which under¬ 
pins viscous flow, therefore providing an unambiguous, 


FIG. 2: Nonlocal response for different resistivity-to- 
viscosity ratios p/rj. Plotted is voltage V(x) at a distance x 
from current leads obtained from Eq. ( |12| ) for the setup shown 
in the inset. The voltage is positive in the ohmic-dominated 
region at large \x\ and negative in the viscosity-dominated 
region closer to the leads (positive values at even smaller \x\ 
reflect the finite contact size a « 0.05u> used in simulation). 
Viscous flow dominates up to fairly large resistivity values, 
resulting in negative response persisting up to values as large 
as p(new) 2 /ri « 120. Nodal points, marked by arrows, are 
sensitive to the p/p value, which provides a way to directly 
measure viscosity (see text). 

almost textbook, diagnostic of the viscous regime. 

Nonlocal electrical response mediated by chargeless 
modes was found recently to be uniquely sensitive to the 
quantities which are not directly accessible in electrical 
transport measurements, in particular spin currents and 
valley currents. |21H23) In a similar manner, the nonlo¬ 
cal response discussed here gives a diagnostic of viscous 
transport, which is more direct and powerful than any 
approaches based on local transport. 

There are several aspects of the electron system in 
graphene that are particularly well suited for studying 
electronic viscosity. First, the momentum-nonconserving 
Umklapp processes are forbidden in two-body collisions 
because of graphene crystal structure and symmetry. 
This ensures the prominence of momentum conservation 
and associated collective transport. Second, while carrier 
scattering is weak away from charge neutrality, it can be 
enhanced by several orders of magnitude by tuning the 
carrier density to the neutrality point. This allows to 
cover the regimes of high and low viscosity, respectively, 
in a single sample. Lastly, the two-dimensional structure 
and atomic thickness makes electronic states in graphene 
fully exposed and amenable to sensitive electric probes. 

To show that the timescales are favorable for the hy- 
drodynamical regime, we will use parameter values esti¬ 
mated for pristine graphene samples which are almost 
defect free, such as free-standing graphene.^ Kine¬ 
matic viscosity can be estimated as the momentum dif- 
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fusion coefficient v ss where 7 ee is the carrier- 

carrier scattering rate, and vf = 10 6 m/s for graphene. 
According to Fermi-liquid theory, this rate behaves as 
7ee ~ (&b T) 2 /Ep in the degenerate limit i.e. away from 
charge neutrality, which leads to large v values. Near 
charge neutrality, however, the rate 7 ee grows and v ap¬ 
proaches the AdS/CFT limit, namely sH/AttIcb where s 
is entropy density. Refs.|121 124] estimate this rate as 
7ee ~ Aa 2 k^T/h, where a is the interaction strength. 
For T = 100 K, assuming Ep = 0 and approximating the 
prefactor as A « ina 123 , this predicts characteristic 
times as short as 7 ” 1 ~ 80 fs. Disorder scattering can 
be estimated from the measured mean free path values 
which reach a few microns at large doping [29]. Using 
the momentum relaxation rate square-root dependence 
on doping, 7 P oc n -1 / 2 , and estimating it near charge 
neutrality, n < 10 lo cm -2 , gives times 7 " 1 ~ 0.5ps, 
which are longer than the values 7 " 1 estimated above. 
The inequality 7 P <C 7 ee justifies our hydrodynamical 
description of transport. 

Momentum transport in the hydrodynamic regime is 
described by continuity equation for momentum density, 

dtPi + djTij = - 7 P pi, Tij = PSij + nViVj + T^\ ( 1 ) 

where T i;1 is the momentum flux tensor, P and y are pres¬ 
sure and mass density, and v is the carrier drift velocity. 
The quantity j p , introduced above, describes electron- 
lattice momentum relaxation due to disorder or phonons, 
which we will assume to be small compared to the carrier 
scattering rate. We can relate pressure to the electro¬ 
chemical potential via P = e J ” ^(n')dn'. Here we work 
at degeneracy, Ep /cbT, ignoring the entropic/thermal 
contributions, and approximating P ss e(n — no)*!*, with 
n the particle number density. While carrier scattering 
is suppressed at degeneracy as compared to its value at 
Ep = 0, here we assume that the carrier-carrier scat¬ 
tering remains faster than the disorder scattering, as re¬ 
quired for the validity of hydrodynamics. Viscosity con¬ 
tributes to the momentum flux tensor through 

= r](diVj + djVi) + (C - rj)dkV k Sij (2) 

where y and £ are the first and second viscosity coef¬ 
ficients. For drift velocities smaller than plasmonic ve¬ 
locities, transport in charged systems is described by an 
incompressible flow with a divergenceless velocity field, 
diVi = 0. In this work, we consider the limit of low 
Reynolds number, yViVj <C y(diVj + djVi), such that the 
role of viscosity is most prominent. At linear order in v, 
we obtain an electronic Navier-Stokes equation 

d t pi - rjV 2 Vi + 7 P pi = -diP. (3) 

This equation describes momentum transport: imparted 
by the external field f = —VP, momentum flows to sys¬ 
tem boundary where dissipation takes place. It is there¬ 
fore important to endow Eq.([3| with suitable boundary 


conditions. In fluid mechanics this is described by the no¬ 
slip boundary condition v = 0. We use a slightly more 
general boundary condition 

nj_ = 0, 7)y = — ac?|| P (4) 

where the subscripts _L and || indicate the velocity and 
derivative components normal and tangential to the 
boundary. The second relation in Eq. © generalizes the 
no-slip condition to account for non-hydrodynamical ef¬ 
fects in the boundary layer on the scales > l = v/"/ ee . 
The model in Eq.Q, equipped with the parameter a, 
provides a convenient way to assess the robustness of our 
predictions. 

It is instructive to consider current flowing down a long 
strip of a finite width. A steady viscous flow features a 
nonuniform profile in the strip cross-section governed by 
the momentum flow to the boundary. Eq.([3]), applied to 
a strip 0 < y < w, yields {—7]d 2 + 7 p mn)v{y) = enE , 
where v(y) and E are the drift velocity and electric field 
directed along the strip (and m is an effective mass de¬ 
fined through the relation p = mnv). Setting a and 
7 p to zero for simplicity, we find a parabolic profile 
v(y) = Ay(w—y), where A = neE /277 and rj = mnv. The 
nonzero shear d v v = A(w — 2y) describes momentum flow 
to the boundary. The net current I = f™ nev(y')dy' = 
(n 2 e 2 /12ri)w 3 E scaling as a cube of the strip width is 
the electronic analog of the Poiseuillc law. Being distinct 
from the linear scaling / oc wE in the ohmic regime, the 
cubic scaling can in principle be used to identify the vis¬ 
cous regime. It is interesting to put the current-field 
relation in a “Drude” form using kinematic viscosity: 
I =» ne Tw wE with t w = w 2 /12v an effective scattering 
time. Evaluating the latter as t w ~ ^(w/vf) 2 7 ee we find 
values that, for realistic system parameters, can greatly 
exceed the naive estimate 7” 1 = w/vf based on the bal¬ 
listic transport picture. This remarkable observation was 
first made by Gurzhi [El- 

Next, we proceed to analyze nonlocal response in 
a strip with transverse current injected and drained 
through a pair of contacts as pictured in FigJT] Un¬ 
like the above case of longitudinal current, here the 
potential profile is not set externally but must be ob¬ 
tained from §. The analysis is facilitated by intro¬ 
ducing a stream function through v = z x Vf, which 
solves the incompressibility condition. At first we will 
completely ignore the ohmic effects, setting a and 7 P 
to zero as above, which leads to a biharmonic equation 
(< d 2 + d 2 ) ip = 0 with the boundary conditions v x = 0 , 
nevy = IS(x) for y = 0 ,w. Using Fourier transform 
in x, we write ip(x,y) = (27t) _1 f dke lkx ipk{y) and then 
determine ipk{y) separately for each k (see Supplemen¬ 
tary Information). Inverting Fourier transform gives the 
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stream function 


ip{x,y) = 


I C dke ikx f e ky + e k( ' w ~ v '> 
ne J 2irik \ e kw + 1 
+afc[ysinlifc(w; — y) + (w — y) sinhfcy]), 


the delta function in the boundary condition for cur¬ 
rent source by a Lorentzian, nev y = Ia/ir(x 2 + a 2 ) at 
y = 0, w. After making appropriate changes in the above 
derivation (namely, plugging e — a l fc l under the integral) 
we find 


where we defined ak = k temh(kw / 2) / (kw + sinhfcw). 
Contours (isolines) of ip give the streamlines for the flow 
shown in Fig[l] While most of them are open lines con¬ 
necting source and drain, some streamlines form loops. 
The latter define vortices occurring on both sides of the 
current path. Numerically we find that vortex centers 
are positioned very close to x = ±w (see Supplementary 
Information). 

We can now explore the electrical potential of the vis¬ 
cous flow. The latter can be found directly from ip(x,y) 
giving 

4>(x, y) = f dke lkx ak[smhk(y — w) + sinhfcy], ( 6 ) 


where we defined (3 = 2rj/(irn 2 e 2 ) (see Supplementary 
Information). As illustrated in Fig[lJ Eq.([ 6 ]) predicts a 
peculiar sign-changing spatial dependence, with two pairs 
of nodal lines crossing at contacts. To understand this be¬ 
havior, we evaluate <p(x, y) explicitly in the regions near 
contacts (x, y) = (0,0), (0, w). Near the first contact, ap¬ 
proximating tanh(fcu;/ 2 ) ss sgnfc, sinhfcy ss le' k ' v sgnk, 
etc, we find 


cp{x, y) « 


PI_ 

2 


dke ikx \k\e~W y 


PI(y 2 - x 2 ) 
( y 2 + x 2 ) 2 


( 7 ) 


(M, \y\ < w). Eq. 0 predicts an inverse-square depen¬ 
dence vs. distance from contacts and also the presence 
of two nodal lines running at ±45° angles relative to the 
nominal current path. Similar behavior is found near the 
other contact, <p(x,y) « We note that 

the r ~ 2 power law dependence is much stronger than the 
lnr dependence expected in the ohmic regime. This, as 
well as multiple sign changes, provides a clear signature 
of a viscous flow. 

The nonlocal voltage measured at a finite distance from 
the current leads (see schematic in Fig{2] inset) can be 
evaluated as V(x) = cp(x,w) — cp(x,0). From Eq.Q we 
predict voltage that is falling off as x~ 2 and is of a neg¬ 
ative sign: 


V(x) 



( 8 ) 


(| a; | < w). Microscopically, negative voltage originates 
from a viscous shear flow which creates vorticity and 
backflow on both sides of the current path, see FigJT} 
Numerically we see that the negative response persists 
to arbitrarily large distances, see p = 0 curve in Fig(2] 
The sign change at very short x, evident in Fig [2] arises 
due to a finite contact size. We model it by replacing 


V(a 


pi 


(x — ia) 2 


cc= -»^) (9) 

(x 2 + a 2 ) 2 [> 


This expression exhibits a sign change at x = a (repre¬ 
senting “the contact edge”) and is negative for all \x\ > a, 
i.e. everywhere outside contacts (this is further discussed 
in Supplementary Information). 

It is interesting to probe to what extent the negative 
response is sensitive to boundary conditions, in particular 
to the no-slip assumption. Extending the above analysis 
to the boundary conditions with nonzero a in Eq.Q we 
find the nonlocal response of the form 


V(x) = PI 


ikx fctanh(fcru/ 2 ) sinhfcui 
kw + (1 + ak 2 ) sinh kw' 


( 10 ) 


where a = 2rja/ne. The expression under the integral 
represents an even function of fc with a zero at fc = 0 
and a symmetric double-peak structure. The peaks roll 
off at |fc| ~ cT 1/2 which sets a UV cutoff for the integral 
similar to that above for a finite-size contact model. Our 
numerical analysis shows that this is indeed the case, with 
a finite a translating into an effective contact size a « 
d 1 / 2 . In other words, the modified boundary conditions 
can alter the response in the proximity to the contact 
while rendering it unaffected at larger distances. 

So far we ignored the bulk momentum relaxation, set¬ 
ting 7 P = 0 in Eq.([3]). We now proceed to show that 
the signatures of viscous flow identified above are ro¬ 
bust in the presence of a background ohmic resistivity 
p = TO 7 p /ne 2 so long as it is not too strong. The dimen¬ 
sionless parameter which governs the respective roles of 
resistivity vs. viscosity is 


e = p(new) 2 /y « 2'y ee 'y p (w/v F ) 2 . ( 11 ) 


For the values 7 ee and quoted above, and taking w = 
1 pm, one obtains e ~ 10. Incorporating finite resistivity 
in the calculation is uneventful, yielding a response 


r dk pe ikx (e qw - l){e kw - l)q 
J irk q + (e qw — e kw ) + g_ [e qw e kw — 1) ’ 


where q 2 = k 2 + ew~ 2 , q± = q ± fc (see Supplementary 
Information). For e = 0 we recover the pure viscous re¬ 
sult, which is negative, whereas for e —► oo Eq. (12 1 gives 
the well-known ohmic result V{x) = £ In [coth( 7 rx/ 2 w)], 
which is positive. With both 77 and p nonzero, the resis¬ 
tance given by Eq.(12) is positive at large x but remains 
negative close to the contact. The dimensionless thresh¬ 
old that determines the possibility of negative electric 
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FIG. 3: Heating patterns for viscous and ohmic flows. 

Viscous flow (panel a) results in a highly complex heating pat¬ 
tern with intense hot spots near contacts and cold arc-shaped 
patches at vortices surrounded by warmer regions. (Also note 
a cold spot at the center where the flow is locally uniform and 
thus W = 0, see Eq. (13)). White arrows show current direc¬ 
tion. Ohmic flow (panel b) shows an essentially featureless 
heat production rate decaying monotonously away from con¬ 
tacts. 


response depends on the actual contact size. For the val¬ 
ues used in Fig. [2] the negative response occurs even at 
fairly high resistivity values corresponding to e < 100 . 

The robustness of negative response can be understood 
by noting that viscosity is the coefficient of the highest 
derivative term in E q.(|3|) and thus dominates at short 
distances x < x* = \Jr\/p(en) 2 . The pervasive character 
of negative response, manifest in Fig(2] will facilitate ex¬ 
perimental detection of viscous transport. The positions 
of the nodes, marked by arrows in Fig[2j vary with the 
ratio p/t 7 , which provides a convenient way to directly 
measure the electronic viscosity. 

The hydrodynamic transport regime also features in¬ 
teresting thermal effects. At a leading order in the flow 
velocity those are dominated by convective heat transfer, 
described by a proportionality relation between entropy 
flux and flow velocity v (to be discussed elsewhere). At 
higher order in v, besides the conventional Joule heat¬ 
ing, vorticity of a viscous flow manifests itself in heat 
production 

W = ?? (diVj + dji’i) 2 = 2ij\(d x +id y ) 2 il>(x,y)\ 2 (13) 

The vorticity-induced heating pattern, shown in Fig(3j 
features hot spots near contacts and cold arc-shaped elon¬ 
gated patches in vorticity regions. In contrast, for an 
ohmic flow the pattern is essentially featureless. This 
makes viscous flows an interesting system to explore with 
the nanoscale temperature scanning techniques. [301 

We finally note that there are tantalizing parallels- 
both conceptual and quantitative-between electronic vis¬ 
cous flows and current microfluidics systems. Namely, a 


model essentially identical to Eq.([3]) describes the low- 
Reynolds (microfluidic) flow between two plates sepa¬ 
rated by the distance h , where 7 p mn = 12 r)/h 2 (it also 
describes viscous electron flow in a 3D slab, see Sup¬ 
plementary Information). A new research area at the 
frontier of nanoscience and fluid mechanics, microfluidics 
aims to manipulate and control fluids at a nanoscale 
with the ultimate goal of developing new lab-on-a-chip 
microtechnologies. Graphene, which can be easily pat¬ 
terned into any shapes without compromising its excel¬ 
lent qualities, can become a basis of electronic microflu¬ 
idics, with multiple applications in information process¬ 
ing and nanoscale charge and energy transport that re¬ 
main to be explored. 
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SUPPLEMENTARY INFORMATION 


Methods 


Here we outline the hydrodynamic approach used be¬ 
low to model electron fluids in the ohmic and viscous 
transport regimes. These regimes differ in the character 
of the constitutive current-field relation and are described 
by different equations for the stream function, harmonic 
in the ohmic case and bi-harmonic in the viscous case. 

For ohmic transport the current-field relation j = crE 
is local, and as a result current is a potential vector field 
with zero vorticity. Indeed the relation E = — V(/>, where 
t f> is electrostatic potential, yields V x j = 0. Relating 
current density and flow velocity, j = nev, and assum¬ 
ing constant particle number density n, we see that the 
velocity field itself is potential. Combining this with the 
continuity equation we write the incompressibility condi¬ 
tion as Laplace’s equation for the electric potential 

-V^ = V 2 0 = O, V 2 =d 2 + d 2 . (14) 

(J u 

Incompressibility allows one to introduce the stream 
function via v = z x V'0 = (—d y ip,d x ip). The stream 
function y) in the ohmic case also satisfies Laplace’s 
equation X7 2 ijj = 0 . 

In contrast, the current-field relation for viscous flows 
is essentially nonlocal. Such flows are described by the 
Navier-Stokes (NS) and continuity equations 

rjS7 2 Vi = neS7 i(j>, = 0, (15) 


where 77 is the dynamic viscosity. Here and below we as¬ 
sume the low-Reynolds limit, as appropriate for the linear 
response regime of interest, and consider linearized NS 
equation. The velocity field of a viscous flow is gener¬ 
ally non-potential since the vorticity ui = V x v = zV 2 t/> 
is non-zero. In this case the stream function is not har¬ 
monic (unlike the Ohmic case) but rather is bi-harmonic, 
satisfying 

(V 2 ) 2 V> = 0. (16) 


We note that the bi-harmonic equation is not conformal 
invariant, in contrast to the harmonic equation in the 
Ohmic case. 

The bi-harmonic equation (16 1 should be endowed with 
the boundary conditions which specify the behavior of 
the flow at system edge. One boundary condition is im¬ 
posed on the normal velocity by the injected and drained 
current I{r ) flowing in and out through the leads. In 
terms of the stream function this reads 




!(r) 

en 


( 17 ) 
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However, since Eq. (16) is fourth-order, a single bound¬ 
ary condition would not suffice and an extra boundary 
condition must be added. Physically, the NS equation 
accounts for momentum exchange among carriers in the 
viscous fluid bulk, whereas the boundary condition must 
account for momentum exchange with the solid bound¬ 
ary. We start with the simplest case of the no-slip bound¬ 
ary condition v; = 0 which we write as 


= o. 


(18) 


Modification for the case of partial slippage is straight¬ 
forward and will be discussed later. With these bound¬ 
ary conditions we seek the spatial dependence of the 
quantities ip, v, <p, and then relate the net current 
/ to the potential difference at remote voltage probes 
V(x) = A(p(x), which defines the nonlocal resistance 
I? n i(:r) = V(x)/I (see inset in Fig.2 of the main text). 


The viscous case 


Here we describe the solution of the problem (151 with 


the no-slip boundary condition in a strip —oo < x < oo, 
0 < y < w. Using the complex variable 2 = x + iy, 


a general solution of the bi-harmonic equation (16) can 


be written in a compact form as ip(x, y) = zf(z) + g(z). 
In the strip geometry the problem can be conveniently 
analyzed in a mixed position-momentum representation. 
Fourier transforming via 


= ^J dke ikx 


V’ k{y ), 


(19) 


the bi-harmonic equation becomes ( k 2 — d 2 ) 2 ipk{y) = 0. 

Here we consider a symmetrical lead arrangement, that 
is the same current profile I{x) at y = 0 ,w. The bound¬ 
ary condition © for the velocity component normal to 
the boundary, gives 

I i,X) 

y')y= 0 ,w — V^)y= 0 ,w — i (20) 

so that 0) = ^k{ w ) = I(k)/enik. The no-slip bound- 
ary condition |l8| ) reads v x (x, 0) = v x (x,w) = 0 giving 
d y ipk{ 0) = d v ipk{w) = 0. Solving for ipk{y) with these 
boundary conditions, and Fourier transforming to posi¬ 
tion space, we obtain 


ip{x,y) = 


dk I(k)‘ 


-.ikx 


+ 


27r neik 

— OO 

ktanhkw/2 


cosh k(y — w/2) 


( 21 ) 


kw + sinh kw 


coshfcic/2 
[y sinh k(w — y) + (w — y) sinh ky\ . 


The flow v = z x V ip, described by this solution, can 
be visualized by the contours ip(x, y) = const which rep¬ 
resent current streamlines. Such streamlines are shown 


by white lines in Fig. la of the main text. Notably, 
there are three clearly distinct groups of streamlines: 
one in the central region, connecting the leads along 
the ‘nominal current path’, and two in the side regions, 
which circle around two points in the middle of the strip, 
r± = (±xo,w/2), where ip(x,y) peaks. These three 
groups of streamlines are separated by two separatrices, 
defined as the outmost streamlines that border the nom¬ 
inal current path region in Fig. la of the main text. The 
two separatrices connect the right and the left edges of 
source and drain, respectively. 

The streamlines circling around the extrema of ip(x, y) 
at the points r± describe vortices generated by the 
flow. The value xq can be obtained from the condition 
d x ip(x, w/2 ) = 0. For the case of point-like current leads, 
I(x) = IS(x), described by a constant I{k) = I , we have 


J dt cos(2uxq/w) ^1 + 


2u tanh u sinh u 
2 u + sinh 2 u 


= 0 , ( 22 ) 


where we defined the integration variable via kw = 2u. 
Solving this equation numerically we find Xq sa w within 
the numerical precision of our method. 


The spatial dependence given by Eq.(21), evaluated for 
point-like leads, I(k) = J, is singular at x —► 0, y —> 0, w. 
The singularity originates from the integrand tending to 
a constant value at large k. The large-A; (“ultraviolet”) 
regularization of the expression in Eq. © can be per¬ 
formed either by making the lead size finite or by alter¬ 
ing the no-slip boundary condition allowing for a small 
partial slip, as discussed below in Section]] Our numeri¬ 
cal analysis indicates that the salient features of the flow, 
pictured in Fig. 1 of the main text, are insensitive to the 
regularization method. 

The electric potential can now be obtained from 


V0(®, y) = ^V 2 [z x ViP{x, y)]. (23) 

ne 

The terms in ip, which are purely exponential in y, are 
harmonic, that is they vanish upon applying V 2 and do 
not contribute the potential cp. The non-harmonic part 
of ip then yields 


<P(x,y) = -^- 2 [ dke ikx 
n(en) 2 J 

x [sinh k(y — w) + sinh ky\. 


fctanh(fc'ic/2) 
kw + sinh kw 


(24) 


The resulting spatial dependence features a remarkable 
behavior which is illustrated in Fig la of the main text. 
The potential <p(x,y), which is zero on the line y = w/2 
by symmetry, exhibits multiple sign changes with several 
nodal lines separating regions where </> is positive and 
negative. For sufficiently small x, i.e. near the nominal 
current path line x = 0, it varies from positive values at 
the source to negative values at the drain. Yet, this de¬ 
pendence is reversed away from the nominal current path. 
In particular, the potential sign everywhere at the strip 













boundaries y = 0 , w outside the leads is opposite to the 
potential at respective leads. This follows from the prop¬ 
erties of the two separatrices of the flow, defined above in 
terms of the borders between the nominal current path 
region and adjacent vortex regions. Consequently, every 
point outside the leads positioned close enough to the up¬ 
per boundary is connected to the symmetrical point near 
the lower boundary by a streamline going opposite to the 
flow in the central region, which translates to a negative 
voltage. This can be also seen by evaluating the voltage 
difference for a pair of points at a distance x away from 
the current leads: 

/ OO 

dke ikx (25) 

-OO 

fctanh(fciy/ 2 ) sinh kw 2q 

kw + sinh kw ’ 7r(en) 2 ’ 


with the parameter [3 defined in the same way as in Eq. ( 6 ) 
of the main text. The voltage V ( x ) is a Fourier transform 
of a positive real-valued function which is even under 
k —> —k. To show that the resulting nonlocal resistance 
R n \{x) = V(x)/I is negative, we first note that the func¬ 
tion under the integral vanishes at k = 0. This means 
that f^° oo V(x)dx vanishes, and therefore V{x ) must be 
negative on a part of the domain —oo<x<oo. Next, 
noting that the quantity under the integral grows as |fc| 
at large fc, 


k tanh(/cw/ 2 ) sinh kw 
k w + sinh kw 


\k\, \k\w » 1 , 


and also taking into account the identity 



1 1 

(x — iO) 2 (x + iO) 2 ’ 


(26) 


(27) 


we conclude that R n i(x) = V(x)/I is negative provided 
x is nonzero and not too large. Indeed, since at \x\ -C w 
the integral (26) is determined by |fc| 1 /w, the identity 
(27) yields 


p , s 2/? 

Rnl(x) «-^ 

X z 


(28) 


For the flow described by Eq. (21) the quantity (d x 
id y ) 2 -il>(x,y) can be evaluated as 


f „ k(y-w/2) 


lk * ) 7_ pup 

\ coshfcw /2 


j,j r 

(d x + id v ) 2 ?l)(x, y) = - / kdke 

ne tt J 

' tanH kw/2 [[1 + k(y - w)}e k y + (1 + ky)e k ^] J . 


kw + sinh kw 


The resulting spatial dependence Q(x,y ), shown in 
Fig. 3a of the main text, is highly nonlocal. The heat¬ 
ing pattern features a non-monotonic dependence with 
multiple cold spots and warmer regions encircling them. 
In particular, there are two cold spots located inside vor¬ 
tices and one positioned at a midpoint between the cur¬ 
rent leads. This is in contrast to the dissipation for an 
ohmic flow which is maximal on the x = 0 line and de¬ 
cays monotonically as x increases (see Fig.3b of the main 
text). This nonlocal behavior of heating can serve as 
another signature of viscous flows. 


The general ohmic-viscous case 

The approach outlined above can be generalized to de¬ 
scribe viscous flows in the presence of ohmic resistivity. 
In this case, momentum transport is governed by 

( d t + 7 p ) Pi - rjS7 2 Vi = -diP, (31) 

where j p describes momentum relaxation due to disorder 
scattering. Using the incompressibility condition and in¬ 
troducing the stream function via v = z x V -0 as above, 
we arrive at 


MV 2 ) 2 - p(en) 2 V 2 ]V’(x, y) = 0, (32) 

where p = 7 p m/ne 2 is resistivity. 

For a strip geometry, we use the representation (19) to 
write Eq.(32) as 


{dy-k 2 ){dl~q 2 )ip k {y) = °, q 2 = k 2 + p(en) 2 /p. (33) 


provided \x\ < w. Lastly, calculating the integral in (26) 
numerically we find that, as x increases, V(x) mono¬ 
tonically grows, remaining negative in the entire domain 
0 < |cc| <00 and approaching zero without changing sign. 
The physical reason for the negative nonlocal voltage is 
a viscous vortex backflow (see Fig. la and accompanying 
discussion in the main text). 

We also quote the result for the rate of heat production 
due to viscous friction: 


Q(x,y) = r] ^2 (diVj+djVi) 2 (29) 

i,j=x,y 

= 2 d[^ly + (i’xx - ^yy) 2 ] = 2.q\{d x + t^) 2 ^! 2 - 


The nonzero p lifts the degeneracy of the eigenvalues, 
leading to a solution which is a sum of four exponents: 


My) = -^ukT, [a± exp(±fcy) + b± exp(±gy)]. (34) 


Writing the boundary conditions (18) and as 

a+e kw + a-e~ kw + b+e qw + b_e~ qw = 1, (35) 

+ fl— + 6+ + b_ = 1, 
k(a+ — a_) + q(b+ — b_) = 0, 
ka+e kw - ka-e~ kw + qb+e qw - 


qb_e~ qw = 0 . 
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and solving for a± and b±, we find 

(e qw - 1) q 


a+ = 


b+ = 


b _ = 


(fc - g) (l - e( fc +«) u ') + (fc + q ) (e«™ - e kw ) ’ 
e fcu ' (e 9 ™ - l)g 

(A; - g)(l - e( fc +«) u ') + (fc + q)(e qw - e kw ) ’ 
(e*™ - 1) fc 

(g -k)( 1 - e ( fc +« , ) w ) + (fc + q)(e kw - e qw ) ’ 
(e to - 1) ke qw 

{q - fc)( 1 - e (fc+?)«') + (fc + g)(e fc, ° - e 9 ™) ‘ 


The potential Fourier harmonic can now be found from 
the relation 


V</> = (—pen + —V 2 ^ v, 
V en ) 


(36) 


Plugging the general solution from Eq.(34) we see that 


only the terms exp(±fcy) contribute, giving 


«*,!/) 


f 


ikx 


a + (k)e ky — a-(k)e~ ky —. (37) 


dk 


k 


Potential difference between the edges of the strip V(x) = 
(j>{w, x) — </>( 0, x) due to ([37| is then evaluated as 


V(x) = 


dk 


e ikx ( e iw 


I ( e kw - 1 )q 


Ip 

7T k q + ( e qw - e kw ) + q- (e^ k + q ) w - l)’ 


{38) 


where q± = q±k. This expression was used to produce 
Fig.2 of the main text. In doing so, we defined the di¬ 
mensionless parameter 


/ \2 P 

e = (enw) — 

1 


(39) 


that characterizes the relative strength of the viscosity 
and resistivity effects, wherein the limiting values e = oo 
and e = 0 correspond to the pure ohmic and viscous 
regimes, respectively. 

As a sanity check, we verify that in these two cases our 
results are in agreement with those found elsewhere. In 
particular, the limit p — >• 0 (i.e. e — > 0 ) can be analyzed 
by expanding the expression under the integral in Eq.( 38) 
in a small difference q — k = \Jk 2 + e/w 2 —q ~ 0(e). Set¬ 
ting e = 0 we find the resistance R n i(x) = V(x)/I which 


is given by (26) and is everywhere negative. Similarly, in 
the opposite limit, 77 —>■ 0 , setting e —> 00 gives qw 00 , 
leading to a± = (e ±kw + l) -1 . Plugging this in Eq.(37) 
leads to an expression 

Ip 


0(x,y) = 


2tt 


f°° c ikx sinh (Hy ~ w/2)) dk 
.oo cosh(fcw/ 2 ) fc 


This gives the nonlocal resistance 
<j>(x, w) — <f>(x, 0) 


^nl(*^) 


(41) 


P 

7T 


r,ikx 


-p— tanh(fcw/2) = — In I coth( 7 rx/ 2 w;)| 
fc 7r 


which is everywhere positive and matches the result oth¬ 
erwise known for the ohmic regime. 

For the general case, with both 77 and p nonzero, we 
expect the resistance I? n i(x) to be positive at large x and 
negative at small x. This behavior can be understood by 


putting Eq. (38) in the form 


V(x) = ^ f°° ? e ikx f(k ), 


(42) 


/(*) = 


a kw 


- 1 


e kw + 1 — | coth (qw)(e kw — 1 ) 


For small enough wavenumbers |fc| <C \fe/w we have 
q = y/pT e/w 2 |fc|. In this case, since |fc|/g -C 1, 
the last term in the denominator of the expression for 
/(fc) can be dropped, giving /(fc) « tanh(fcw;/ 2 ) and 
producing an expression for V(x) which is identical to 
Eq. (41). This means that at distances larger than ^ 
the quantity R n \(x) exhibits the ohmic behavior, Eq.(411, 
and is therefore positive. At the same time, for large 
enough wavenumbers \k\w 1 we have qw 1 and 
coth (qw) « 1. In this case we can simplify the expres¬ 
sion for /(fc) to read 


f(k) 


sgn fc 


q(q+ \k\)w 2 . 


(43) 


This expression can be used to describe the behavior 
of R nl (x) at | a" | < w. The latter is particularly sim¬ 
ple for e <C 1 (i.e. at low resistivity or high viscos¬ 
ity). In this case |fc|w -C 1 implies q ~ |fc|, giv¬ 
ing /(fc) ss sgn (k)k 2 w 2 /e. Plugging this expression in 
Eq. ( [42] ) we obtain the result identical to that in the pure 
viscous case, i.e. negative R n \(x) = —2/3/x 2 given in 
Hq. |28j) above. A more complex behavior is found for 
e > 1 (corresponding to high resistivity or low viscos¬ 
ity), giving /(fc) that behaves differently in the domains 
1 <C \k\w -C y/e and |fc|tu yfe. Namely, 


m = 


sgn fc 


1 < |fc|w < yfe 


s -^k 2 w 2 \k\w > yf€ 


(44) 


In the first case, after plugging in Eq. (42), we find the 


logarithmic dependence V(x) = (2/-K)pI\n(w/\x\), giv¬ 
ing R n \(x) of a positive sign. In the second case, we find 
the familiar viscous spatial dependence R n \(x) = —2/3/x 2 
of a negative sign. The two behaviors are restricted to 
the domains 7 )= < |x| < w and |ar| < 7 )=, respectively. 

The above represents a fairly dramatic behavior, which 
is simultaneously non-monotonic and sign-changing: as 
x approaches the leads R n \(x) first grows to higher 
and higher positive values, behaving as R n i(x) = 
(2/Tt)p\a.(w/\x\), and then abruptly drops and reverses 
its sign, behaving as R n i(x) = —2/3/x 2 . The sign change 
takes place at |x| « The complexity of this spatial 
dependence can be linked to the fact that viscosity repre¬ 
sents a singular perturbation of transport equations (as 
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a coefficient in front of the highest derivative). As a re¬ 
sult, even when viscosity is small, it always dominates at 
sufficiently short distances. 

Further, one can argue that the sign change of i? n i(x) 
occurs at |x| ~ ^ also when e <C 1 (i.e. at low resistiv¬ 
ity or high viscosity). Indeed, the behavior at distances 
|x| < ^ derives from wavenumbers \k\w y/e. In this 
case, q ss k for positive k and q ss — k for negative k. 
Assuming without loss of generality k > 0 and expand¬ 
ing Eq. (381 in small q — k (while allowing the quantity 
kw to be either large or small) we obtain an expression 
identical to Eq.( 26) found in the pure viscous case. This 
expression was investigated analytically and numerically 
in the main text and shown to produce R n i(x) of a neg¬ 
ative sign. 


The robustness of the negative nonlocal resistance 


It is instructive to compare the behavior A. n i(x) = 
—2/3/x 2 , which was obtained above assuming point-like 
current leads and no-slip boundary conditions, to that 
found under more general assumptions. One interesting 
generalization involves altering the boundary conditions 
to allow partial slippage at the boundary with the veloc¬ 
ity proportional to the electric field: 


v x (x, y ) v =o = -d y 4>(x, y) v = 0 = a(x)d x <j>(x, y) y =o, (45) 


and similar for y = w. Here for the sake of generality 
we made the surface slip factor a(x) position dependent, 
which can serve as a model of structural or chemical mod¬ 
ulation at system boundary. For a constant a , the solu¬ 
tion can be obtained by generalizing the approach out¬ 
lined above. This is done by seeking the stream function 
in the form |l9| ). Potential has the same y-dependence as 
above, 4>k{y) = 6*[sinhfc(y—to)+sinhfcy], where the coef¬ 
ficients b k must be found from (23) and the new boundary 


condition. After some algebra we find 


bk = 


iv 


fctanh(fcw/2) 


7r (en) 2 kw + sinh/cui + 2(rja/ne)k 2 sinh kw 


(46; 


This expression was used in deriving Eq.(10) of the main 
text. The a-dependent term in the denominator changes 
the large -k asymptotic of bk- This translates to the x 
dependence in which the negative singularity — 1/x 2 is 
replaced, at small enough x, by a much weaker positive 
singularity ln(l/x). The value and sign of the nonlo¬ 
cal resistance J£ n i(x) remain unaffected (and negative) 
at not-too-small distances x yj rja/ne , indicating that 
the effect of nonzero a is inessential and can be neglected. 
We note parenthetically that it may be interesting to ex¬ 
ploit position dependence a(x) by adding a small periodic 
modulation via a(x) = a + 5a cos(fcx). This may lead to 
a new behavior characterized by a modulation </(x) with 
the same periodicity, i.e. a chain of vortices. 


In a realistic setting, the singular behavior I? n i(x) ~ 
— 1/x 2 is also regularized by a finite contact size. There 
are two main types of contacts used in the measurements 
on high-mobility graphene: 

1) ideal metal contacts; 

2) narrow graphene channels shaped through etching 
so that they connect seamlessly to the graphene bulk. 


The effect of a spread-out current be modeled by tak¬ 
ing the current I(x) to be distributed over a finite cross- 
sectional area. In this case Eq. (26) reads 


/ GO 

dkl(k)( 

-GO 


ikx k tanh(fcu>/2) sinh kw 
kw + sinh kw 


(47) 


For currents spread over a region of size £ we expect the 
potential at |x| i to remain unaffected, approximately 
following the —1/x 2 dependence as |x| decreases. How¬ 
ever, as |x| approaches £, we expect the potential to re¬ 
verse its sign and become positive at the contact. As a 
crude model, we illustrate this effect in the main text us¬ 
ing a Lorentzian distribution, in Fourier representation 
described by I(k) = I exp(—a|fc|). 

A more realistic model should account for the finite 
size and sharp edges of the contacts. Here we discuss the 
case of small-size contacts such that £ -C w, a practically 
relevant case which is fairly straightforward to analyze. 
For x <C w, using the \k\w 1 asymptotic given in 


Eq.(26) and the identity (271, we can re-write Eq.(47) in 


the position space as follows: 

/ OO 

dk I (k)e ikx \k\ 

-OO 


(48) 


' —OO 

poo 


= -13 


dx' 


I(x') 


(x — x' + i0)- 


+ c.c. 


At first, with the model 2 in mind, we ignore the equipo- 
tential condition and study potential obtained from a 
fixed current distribution. The resulting behavior can 
be exemplified by a current density constant within the 
region — i < x < £, representing contact. In this case 


Eq.(49) predicts voltage 


V(x) = - 


2 13 


.2 _ (2 


(49) 


which is negative at |x| > £ and reverses sign at the edges 
x = ±A the second-order pole x~ 2 found above for point 
contacts is now split into two first-order poles at x = ±£. 
The large-x asymptotic l/(x 2 — £ 2 ) = 1/x 2 + 0(£ 2 /x 4 ) 
then yields the dependence far outside contacts, 

V(£ < \x\ < w) = —^7, (50) 

x z 

that matches the expression found for point-like leads. 
Deviation from this behavior occurs in the small space 
region x ~ £, as expected. 
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One can further smoothen the singularity in V ( x ) by 
making the current vanish at the edge. For example, this 
is the case for the parabolic distribution 


0,w 


3I(£ 2 -x 2 )/4£ 3 , 


0, 


1*1 < 

lari > 


(51) 


After Fourier transforming we find I{k) = 6/[sin(/c£) — 
klcos(k£)\. Since I(k) tends to a constant / when 
kl —> 0, the behavior at x £ remains unaffected by 
the details of current distribution at the leads. In the 


case (511, integral (49) is conveniently evaluated by writ¬ 
ing {x _J +i0) 2 = and moving the derivative 

on /( x') via by parts integration: 


V{x) = p f dx' f 

J —OO V 


= W [ dx' 


I'{x') 


+ c.c. 


2£ 3 


x — x' + iO 

/ 

+ C.C. 


(52) 


x 


x — x' + iO 


Integrating over x' with the help of the indefinite integral 
f ^ U a^u ~ f ^ u U ~a-v, a = ~ u ~ a l n ( a — u ) gives a closed- 
form expression for V (ar) which is valid both outside and 
inside the leads: 


£ 3 V \x + £\ 


(53) 


The behavior at large x , obtained from the asymptotic 
xln = — 2£ — + 0(£ 5 /ar 4 ), agrees with the depen¬ 


dence far outside the leads found above, see Eq. (501. 


Notably, the voltages given in Eq.(49) and Eq.(53) are 


not constant inside the interval — l < x < £. This behav¬ 
ior, which is physical for model 2, does not describe ideal 
metal contacts (model 1). In the latter case we expect 
the potential to be constant within the leads. To satisfy 
the equipotential condition one has to find the current 
distribution I(x) and potential V{x) self-consistently, in 
a way that ensures that the resulting V (x) does not vary 
inside the leads. This can be accomplished by treating 
the relation between the potential and current, 


V(x) =pj dx' ^ 


I'(x') 


x — x' + iO 


+ c.c. 


(54) 


The electrostatic potential spatial dependence, which in¬ 
cludes the contributions due to the external field Eo and 
the charges it induces on the strip, can be written as 


y) = - Re A(C 2 - £ 2 ) 1/2 , 


(56) 


where we introduced a complex variable ( = x + iy. This 
expression describes a function that takes a zero value 
within the strip —£ < x < £ and has an asymptotic 
behavior of the form <f>(x,y) ~ — Xx at distances large 
compared to £. The charge density on the strip can be 


found from Eq. (56) with the help of Gauss’ law, giving 


r{x) = 


Xx 


1, 


|a;| < £ 

2tt(£ 2 - a; 2 ) 1 / 2 ] 0, \x\> £ 


(57) 


Potential due to a(x) accounts for the difference between 
the uniform field potential 4>o(x,y) = —Xx and the po¬ 


tential given in Eq. (56). For points in the y = 0 plane 
this gives 

/ OO 1 

cr(a; / )21n -- -dx' = Xx— Re X(x 2 — i 2 ) 1 ^ 2 , (58) 

-oo \x-x'\ 

where the branch of ( x 2 — l 2 ) 1 / 2 at negative x is found 
by analytic continuation from positive x, giving ( x 2 — 
l 2 y/ 2 S gnx. Taking the derivative with respect to x we 
find the electric field x component E x (x) y ~ oi which gives 
the relation 


t(x')- 


X — X' 


-dx' = —X + Re 


A|a 


(a; 2 — t 2 ) 1 / 2 ’ 


(59) 


Comparing this to the integral equation (|54j) we can re¬ 
late current density in the leads and the potential V(x). 


In doing so we note that the right-hand side of Eq. (59) 


is constant for — £ < x < £ since the last term in (59) 
vanishes for such x. We therefore see that a solution 


of Eq. (54) can be obtained by identifying a{x) with the 
quantity —pi'(x), and A with the potential Vo value 
within the lead. 

Once the correspondence between the two problems is 
established the distribution of current in the lead can be 
found by solving the equation 


as an integral equation for the function I{x) localized in 
the interval — t < x < £. 

Solution of this integral equation under the condition 
that V (x) takes a constant value Vo f° r all — £ < x < £ can 
be found by making use of an auxiliary electrostatic prob¬ 
lem, chosen so that the potential V (x) and the derivative 
of the current I'(x) translate into the electrostatic field 
and charge density, respectively. To that end, we con¬ 
sider an ideal conducting strip of width 2£ in 3D placed 
in a uniform external electric field Eo = Ax. The strip is 
taken to be infinite-length and zero-thickness, and posi¬ 
tioned in the y = 0 plane, such that 

— £ < x < £, y = 0, —oo < z < oo. (55) 


Pl'(x) = —a(x). (60) 


Integrating over x gives the current density which is of 
a constant sign within the lead and vanishes at the lead 
edges: 

nx) = ^(£ 2 -x 2 ) 112 , -£<x<£. (61) 

Evaluating the net current as Iq = J^ e I(x)dx = 
P~ 1 Xtt£ 2 /2 and setting A = Vo, we obtain the contact 
resistance of the lead 


Rc = 


Vo 

Io 


4p 

£ 2 ' 


(62) 
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The inverse-square relation R c oc £~ 2 is sharply distinct 
from the log dependence R c oc In (w/£) found in the ohmic 
regime, and can therefore serve as a hallmark of the vis¬ 
cous regime. The potential outside the leads can now be 
found simply as the right-hand side of Eq. (59): 


V(\x\ >£) = V 0 ~ — 
\x 

_ 4 P_ 

( 2 . 2 _£ 2 ) 1 / 2 (( ; e 2 _ 


VqM 

-* 2 ) 1/2 


£ 2 y/ 2 + \x\) 


lo¬ 


rn 


The large-a; asymptotic matches the dependence V(x) = 
—2/3/o/x 2 , as expected. Interestingly, V(x) takes a con¬ 
stant positive value Vq within the source lead, while being 
negative outside the lead and exhibiting a square-root di¬ 
vergence as x approaches the edges x = ±£. The behavior 
on the drain side is opposite to that on the source side. 
The origin of the negative divergence can be understood 
by considering the streamlines outside of the nominal cur¬ 
rent path region, which represent closed loops (see Fig. la 
of the main text). For example, a streamline approach¬ 
ing one of the edges of the source lead along the strip 
boundary y = w must turn and leave along the respec¬ 
tive separatrix which borders the nominal current path 
region. Since the separatrices make finite angles with 
the strip boundaries, the closer the streamline is to the 
boundary, the closer it comes to the lead edge and the 
sharper the turn it must take. Obviously, to take a sharp 
turn, a large enough electric field is required. This sin¬ 
gular behavior may change upon introducing partial-slip 
boundary conditions, to be discussed elsewhere. 


To derive the above result, we consider viscous flow of 
an electron fluid in an infinite slab 


— oo < cc < oo, 0 < y < w, 0 < z < h, (64) 


in which current is injected and drained through the leads 
placed at the slab edges y = 0,w. We assume for sim¬ 
plicity that the current density at the leads has a ver¬ 
tical parabolic profile proportional to z(h — z). In this 
case, the flow inside the slab is planar, v = (v x ,v y ), 
taking on a globally uniform ^-dependence: Vi(x,y,z ) = 
Ui(x,y)z(h — z)/h 2 . The assumption of a parabolic z- 
dependence does not impact the generality of our anal¬ 
ysis since a generic profile at the leads will transform to 
the parabolic profile in the vicinity of the leads, resulting 
in a flow in the system bulk with a parabolic z depen¬ 
dence at the distances from the leads exceeding h. The 
parabolic model is therefore expected to yield fully accu¬ 
rate results for the slab widths w h independent of the 
contact size, providing also a reasonably good approxi¬ 
mation for w h. 

To transform the 3D linearized stationary NS equation, 


[q(dl + d 2 + d 2 ) - p(en) 2 ]vi = neV i<j>, (65) 


to a 2D NS equation, we integrate Eq. (65) over 2 taking 


into account that the potential <f> is ^-independent be¬ 
cause the flow is planar. Integrating the parabolic profile 
of the flow velocity over the cross-section then yields 


[q{d 2 + d 2 ) - p{en) 2 - 12 y/h 2 ]ui = 6 neVj 0 . ( 66 ) 


Three-dimensional viscous electron flow 

Here we discuss the possibility to observe current vor¬ 
tices and negative nonlocal response in viscous charge 
transport in a 3D conducting slab of a small but finite 
thickness. In this case, as we will see, our transport equa¬ 
tions are isomorphic to those describing the low-Reynolds 
(microfluidic) flow between two plates separated by a 
distance h, where the effective resistivity is defined by 
p(en) 2 = 12 p/h 2 . Comparing this to the dimensionless 
control parameter introduced above, e = p{new) 2 / 77 , that 
governs the respective roles of resistivity versus viscosity 
in a strip, we see that for microfluidic flows between two 
plates of width w, separated by a distance h , this param¬ 
eter equals e = 12w 2 /h 2 . 


Taking the two-dimensional curl of (661 we obtain trans¬ 
port equations identical to those for the mixed ohmic- 
viscous model, Eq. (32), albeit with a renormalized resis¬ 
tivity: 


p(en) 2 —» p(en) 2 + 12 p/h 2 . (67) 

It is instructive to compare this with the dimensionless 
threshold that determines the possibility of negative elec¬ 
tric response. For the contact size parameter used to 
produce Fig.2 of the main text, we found the threshold 
value e = p(enw) 2 /r] ss 120. For a slab of thickness h 
this translates into p(enw) 2 /p + 12w 2 /h 2 < 120. This 
condition is somewhat more strict than that in 2D lay¬ 
ers, however it can still be met even when the two slab 
dimensions are equal, w ~ h, e.g. in a wire. 








